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Abstract 



We formulate a Swendsen- Wang-like version of the geometric cluster algo- 
rithm. As an application, we study the hard-core lattice gas on the trian- 
gular lattice with the first- and the second-neighbor exclusions. The data 
are analyzed by finite-size scaling, but the possible existence of logarithmic 
corrections is not considered due to the limited data. We determine the crit- 
ical chemical potential as (i c = 1.75682(2) and the critical particle density 
as p c = 0.180(4). The thermal and magnetic exponents yt = 1.51(1) ~ 3/2 
and yh = 1.8748(8) « 15/8, estimated from Binder ratio Q and susceptibility 
X, strongly support the general belief that the model is in the 4-state Potts 
universality class. On the other hand, the analyses of energy-like quantities 
yield the thermal exponent yt ranging from 1.440(5) to 1.470(5). These values 
differ significantly from the expected value 3/2, and thus imply the existence 
of logarithmic corrections. 
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I. INTRODUCTION 



Lattice gases, together with the Potts (including Ising) and the O(n) model, play an 
important role in the statistical mechanics. They are used to describe universal properties 
of many complex physical systems, ranging from simple fluids to structural glasses and 
granular materials. Lattice-gas models are generally defined as follows. For a given lattice, 
a number of particles is randomly distributed over its vertices with the constraint that each 
vertex can at most be occupied by one particle; the density of particles is controlled by the 
chemical potential //. Particles on different vertices can interact with one another-normally 
through two-body interactions. Accordingly, the reduced Hamiltonian (already divided by 
kT with Boltzmann factor k and temperature T) of a lattice-gas model can be written as 



where a = 0, 1 represents the absence and the presence of a particle, respectively. The second 
term with amplitude Kin describes the first-neighbor interactions, and the third one with 
K 2N is for the second-neighbor couplings; further-neighbor interactions can be included, as 
denoted by symbol • • •. 

In the study of lattice gases, one often takes the hard-core limit: in Eq. (1) the couplings 
K outside a certain range (r > ro) are set zero, while Kir < ro) is taken to the limit 
K — > — oo; namely, the particles have a hard-core of radius r . A particular example is the 
lattice-gas model with nearest-neighbor exclusion: K 1N — > — oo while all further-neighbor 
couplings are zero. Unlike the Potts model and the O(n) spin model, the nature and the 
universality of the phase transitions in lattice-gas systems depend on the lattice structures. 
For instance, the lattice-gas model with nearest-neighbor exclusions on the square and the 
honeycomb lattice is believed to be Ising-like, while that on the triangular lattice (Baxter's 
hard- hexagon model [1,2]) belongs to the 3-state Potts universality class. 

Extensive investigations have been carried out for lattice gases, and many theoretical 
and numerical approaches are applied. This includes exact calculations (mainly by Baxter 
and coauthors), series expansions (like high-temperature and low-temperature expansions), 
cluster variation method, transfer matrix calculations, and Monte Carlo simulations etc. The 
critical free energy of Baxter's hard- hexagon lattice gas was exactly calculated [1,3-5]; the 
critical chemical potential is known as p c = ln[(ll + y/5)/2], and the critical particle density 
is p c = (5 + \/5)/10. Baxter's hard-square model [1,3-5], defined by Eq. (1) on the square 
lattice with K^ — > — oo and finite K 2 n, is known to have a tricritical point (/x tc , K 2 N,tc) 
in the tricritical Ising universality class; the tricritical point lies at p tc = — ln[8(l + a/5)], 
K 2 N,tc = In (3 + y/5i), with p tc = (5 + \/5)/10. Using the transfer-matrix technique, Guo 
and coauthors [6] determine the critical point of the hard-core square lattice gas up to the 
eleventh decimal place, /i c = 1.334 015100 277 74(1), p c = 0.367 742 999 0410(3). Recently, 
Monte Carlo simulations were carried out for square lattice gases with the hard-core radius 
up to the fifth neighbors [7]. The nature of phase transitions was found to be continuous 
for exclusions up to IN, to 2N, and to 4N, and to be discontinuous for exclusions up to 3N 
and 5N, where symbols iN represents the ith neighbors. 

In this work, we shall consider the hard-core lattice gas on the triangular lattice with 
the first- and the second-neighbor repulsions. The triangular lattice in this case can be 
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divided into four sublattices (see Fig.l), and for sufficiently high density particles prefer to 
occupying one of the four sublattices. Thus, one would expect that, if it is second order, 
the melting of ordered phase should be in the 4-state Potts universality class. However, 
several studies at the end of 60s last century suggested that the phase transition is first 
order [8,9]. Later, Bartelt and Einstein [10] reexamined this model by a phenomenological 
renormalization-transfer- matrix scaling. The largest system size in their study is 16 x oo. 
Slowly convergent finite-size corrections were observed. It was estimated that the thermal 
and the magnetic critical exponents are y t = 1.400 and y h = 1.885, respectively. Despite the 
noticeable deviations from the exact values y t = 3/2 and y h = 15/8 [11,12], these estimates 
are in favor of the 4-state Potts universality in view of the possible occurrence of logarithmic 
corrections. 

Here, we aim to provide an independent study of this model by means of Monte Carlo 
simulations. This seems justified since no rigorous argument exists about the nature of 
the phase transition and the evidence in Ref. [10] is not very strong. To properly analyze 
finite-size corrections, more accurate numerical data, particularly for large system sizes, 
are desirable 1 . Our task becomes now feasible because of the availability of efficient cluster 
algorithm for lattice-gas models-the geometric cluster algorithm-and the rapid development 
of computer industry in the past few decades. The geometric cluster algorithm [15,16] moves 
round a fraction of particles over the lattice according to geometric symmetries, such as the 
spatial inversion or rotation symmetries of the triangular lattice; detailed description will 
be given in Sec. II. The algorithm does not change the total number of particles, and it is 
combined with the Metropolis steps in order to simulate lattice-gas systems in the grand- 
canonical ensemble. In comparison with simulations using the Metropolis method only, 
critical slowing down is significantly suppressed. Therefore, we are able to simulate systems 
as large as 400 x 400 within reasonable computer resources. 

II. GEOMETRIC CLUSTER ALGORITHM AND SAMPLED QUANTITIES 

A. Geometric cluster algorithm 

The geometric cluster algorithm was first proposed by Dress and Krauth [17] in the study 
of hard-core gases in continuous space. Unlike the well-known Swendsen-Wang (SW) cluster 
method which flips spins, the elementary operation in this algorithm is to move particles. For 
hard disks, the percolation threshold of the cluster formation process deviates significantly 
from the phase transition of the model. This is unfortunate since it affects the efficiency of 
the algorithm. 



1 Sometimes, when finite-site corrections are not properly taken into account, wrong conclusions 
can be reached. For instance, from the Metropolis simulations of the lattice-gas model on the 
simple-cubic lattice with the first-neighbor repulsions, Yamagata estimated the critical exponents 
as (3/j = 0.311(8) and 7/1/ = 2.38(2) [13], which would imply y h = 2.689(8). This result is 
significantly different from the general accepted value yh = 2.4816(2) for the Ising universality 
class in three dimensions [14]. 
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A single-cluster version of the geometric cluster algorithm was later developed by Heringa 
and Blote [15,16] for lattice models like the Potts model and the lattice gases. Here, we shall 
briefly describe it in terms of the lattice-gas model (1) with a soft-core radius of a lattice 
unit (finite K := K 1N < and all other couplings are zero) on the square lattice with 
periodic boundary conditions. For such a geometry, one can set a Cartesian coordinate by 
taking any two perpendicular lines of lattice sites as the x and the y axis, respectively. It 
can be seen that the Hamiltonian of the system is invariant under geometric operations 
like reflections about the x or the y axis or inversion about the center of the coordinate. 
Further, any configuration will be restored if an operation is subsequently applied twice- 
namely, these operations are self-inverse. One can employ any of such geometric operations 
to formulate a cluster algorithm. Let a pair of nearest-neighboring sites i, k be mapped onto 
i', k', respectively. One denotes the energy difference when a neighbor k of % is interchanged 
with k! as A ik , which is A^ = K(a,i<j k + oyov — — 0^0^). The algorithm then involves 
the following steps: 

1. Choose a random site i: both i and i' belong to the cluster. 

2. Interchange Oi and oy. 

3. For all neighbors k of % that do not belong to the cluster yet, do the following: 

• If Ajfc > 0, do the following with probability p — 1 — exp(— A ik ): (a) interchange 
a k and <r k > (k and k' are included in the cluster), (b) write k in a list of addresses 
(called the stack). 

• If A ik < 0, do nothing. 

4. Read an address j from the stack. Substitute j for i, and execute Step 3. 

5. Erase j from the stack. 

6. Repeat Steps 4 and 5 until the stack is empty. 

When the stack is empty, the cluster is completed. Since the elementary operation is to 
interchange spins and oy, the total number of particles does not change in the algorithm. 

For the above geometric cluster steps, the detailed balance has already been proved 
[15,16]. The efficiency of this algorithm for different models has also been demonstrated. 
For the Ising model, it was shown that the percolation of the formed clusters coincides 
with the thermal phase transition, reflecting the efficiency of the algorithm. In fact, in the 
canonical ensemble (the total number of particles is fixed), it can be shown that, for many 
models, no critical slowing down exist for some quantities [18]. 

Here we shall formulate a full-cluster version of the geometric cluster algorithm in an 
analogous way as the Edward-Sokal picture for the well-known SW cluster method for the 
ferromagnetic Potts model [19,20]. We consider the lattice-gas model (1) with finite nearest- 
neighbor interactions [K := < 0) on a one-dimensional chain with sites labelled as 
i = ±1/2, ±3/2, ±5/2, • • •. Instead of writing the Hamiltonian for a fixed number of particles 
as a sum of the nearest-neighbor couplings like in Eq. (1), we rewrite it as 
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H\n v =n = ^2 7~Li = —K E ( ai °' i + 1 + cr-iC-i-i), (2) 

i=l/2 i=\j1 

where N p = Y2i a i an d the constant N denotes the total number of particles. The Hamil- 
tonian (2) is obtained by applying the reflection about the center % = 0-a geometric op- 
eration. In this form, the 'building blocks' of the Hamiltonian is no longer a pair of 
neighboring sites, but two pairs of them. If one only uses the spin-interchange operation 
rj := <7j <-> (T_j, the energy associated with each building block is of two levels at most: 
£i(a) := —K((Jiai +1 + a_ i o"_ i _ 1 ) and E^ip) := —K{piO-i-\ + <7_j<7j + i). The former £ 1 refers 
to the status that no operator r\ is applied or it is applied at both vertices i and j; instead, 
the latter £ 2 means that rj is applied at i (or j) only. The values of £\ and £2 depend on 
the spin configuration a on the building block. For the lattice-gas mode (2), these values 
are shown in Table I. Let us denote the lower and the upper level of £\ and £2 as £\ ow and 
£ upP , respectively. The statistical weight associated with each block in Eq. (2) reads 

e -«i(*) = e -£app(i + Vi S £lAow ) with ( Vi = e e »™- £l ™ - 1). (3) 

On this basis, the partition sum becomes 

00 00 

z w=N = e n e ~ £upp{ff) n (! + ^j. (4) 

{cr}:AT p =7V i=l/2 i=l/2 

Analogously as mapping the Potts model onto the random-cluster model, one introduces a 
bond variable bi to graphically represent the expansion of the second product in Eq. (4): 
if the term Vi5s lt £ lov/ is taken, an occupied bond bi = 1 is placed between sites i and i + 1; 
otherwise, no bond is placed (6, = 0). This leads to a joint model 

00 

Z\n v =n= E II e_£uPp(CTl EMfcAJ*. (5) 

{a}:Afp=Ar i=l/2 {b} 

where the second sum is over all possible bond configurations that are consistent with the 
spin configuration, and we have already assumed the conventional symbol 0° = 1. Given a 
spin configuration {a}, the expression (5) allows us to place bonds and construct clusters 
as follows: if the spin configuration on a block i is at the lower-energy level £\ ow , one places 
a bond h = 1 with probability fj/(l + i>j); otherwise, place no bonds. Note that a bond 
connects four lattice sites, since it is placed on the blocks. Lattice sites connected through 
a chain of occupied bonds form a cluster. The condition 5s 1: £ lovr = 1 for a block can be 
hold either by doing nothing or interchanging both spins ((J, <-> o"_j,o"j + i <-> <r_j_i). Thus, 
for a spin-jointed-bond configuration, for each cluster one has the freedom to choose the 
do-nothing or the spin-interchange operation, and apply it to all lattice sites within the 
cluster. Accordingly, a Swendsen- Wang-like geometric cluster algorithm can be formulated 
as follows. 

1. Choose a geometric transformation such that every building block in the Hamiltonian 
consists of two pairs of neighboring couplings and the associated energy is only of two 
levels under the spin-interchange operation. 
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2. For each building block % (containing four lattice sites), if its spin configuration is at 
the lower-energy £ low , place a bond with probability Pi — Vi/(1 + i>j); otherwise, place 
no bond. 

3. Construct clusters according to the occupied bonds. 

4. Independently for each cluster, randomly choose the do-nothing or the spin-interchange 
operation with probability 1/2; apply the chosen operation to all lattice sites within 
the cluster. 

A Monte Carlo step is completed, and a new spin configuration is obtained. We expect that 
the analogy between our formulation of the geometric cluster algorithm and the well-known 
SW method can help the reader to understand the geometric cluster algorithm. 

We conclude this subsection by mentioning the following. (1), like the SW method, 
the essence in the geometric cluster process is that the energy of each building block has 
two levels only under the spin-interchange operation (the energy of a building block may 
have more than two levels if other operations-like the spin-flip operation-are allowed). (2), 
normally, the rewriting of the Hamiltonian as a sum of proper building blocks is obtained 
by applying some global geometric transformation, such as the inversion about the center 
and reflection etc. (3), in the canonical ensemble, if one uses the geometric cluster algorithm 
only, a large number of spatial transformations should be available such that each lattice 
site in a system should be able to reach any other lattice site by a finite number of geometric 
mappings. For the torus geometry, this can be easily achieved since any site can serve as 
the center of the aforementioned Cartesian coordinate. In case that the geometric cluster 
method is itself non-ergodic, it can be combined with other algorithms like the Kawasaki 
dynamic. (4), for simulations in the grand-canonical ensemble, other Monte Carlo methods 
have to be used. 

B. Sampled quantities 

For the lattice-gas model (1), the triangular lattice is divided into four sublattices. The 
particle density is then sampled as 

p'" = £ £ * («) 

where V = L x L is the volume of the lattice and the sum is over each sublattice, labelled 
as i = 1, 2, 3, 4. The global particle density is then p = (pW + p^ + p< 3 > + p (4) )/4. On this 
basis, we measured the second and the fourth moment of the magnetization density as 

M 2 4tt (P (t) - P {3) ? and M * = i M ") 2 . ( 7 ) 

i=l j=i+l 

where factor 1/3 is for normalization purpose such that Ai 2 is a unity when the chemical 
potential is infinite-one of the four sublattices is fully occupied. The magnetic susceptibility 
is x — V{M. 2 ). Then, we define a dimensionless ratio as 
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O - {Mr (8) 

This ratio at criticality approaches a universal value for L — > oo, and is known to be very 
useful in estimating critical points. 

Since a pair of first- (or second-) neighboring sites cannot be both occupied by particles, 
we sampled the third- neighbor correlation as an energy-like quantity 



^7 E °^ ( 9 ) 



6V 

{ij}e3N 

where the sum is over all the third-neighbor pairs. Correspondingly, a specific-heat-like 
quantity is defined as C e = V((S 2 ) — (£) 2 )- We also measured the compressibility C v = 

vm-( P n 



III. RESULTS 

Using a combination of the Metropolis and the geometric cluster algorithm, we simulated 
the lattice-gas model on the triangular lattice with first- and second-neighbor repulsions. 
Periodic boundary conditions were used in the rhombus geometry shown in Fig. 1. System 
sizes took 15 values in range 8 < L < 400. Several geometric cluster steps are performed 
between subsequent Metropolis steps, such that the total number of particles moved by the 
former is also approximately equal to V. Significant critical slowing down was observed, 
with a dynamic exponent about z ~ 1. 

Some primary simulations for relatively small system sizes were first carried out to find 
the approximate location of the critical point from the intersection of the Q data for different 
system sizes L (we were also guided by the result fi c ~ 1.7599 in Ref. [10]). Then extensive 
simulations for large sizes were performed near /z c = 1.757. 

Figure 2 shows parts of the Monte Carlo data of the dimensionless ratio Q. According 
to the least-squares criterion, the Q data were fitted by 

Q(/i, L) = Q C + qi (n - ii c )L yt + q 2 (/2 - fi c )L 2yt 

+6iL Wl + b 2 U» + b 3 L y3 + ri L y r + Cl (pL- pi c )L yt+ yi , (10) 

where g« , 6j , q , and are unknown parameters, and Q c is the universal value. The 
terms with qi describe the contributions of the thermal field due to deviation from the 
critical point, those with exponent yi account for finite-size corrections, and the one with 
Ci is for the mixed effect of the relevant and irrelevant thermal fields. The term with 
exponent y r = d — 2yh arises from the regular part of the free energy, where the magnetic 
exponent yn was fixed at 15/8 for the four-state Potts model. The detailed derivation of the 
finite-size scaling formula (10) can be found in Ref. [14]. In principle, one should include 
logarithmic corrections in Eq.(10) [21], since the transition is expected to be in the 4-state 
Potts universality class. Unfortunately, the limited system size does not allow us to include 
such corrections (accordingly, the statistical error margins of our following results should 
be taken carefully). From the numerical results for the tricritical 4-state Potts model [18] 
where the marginal field is absent, we learn that there exist correction terms with exponent 
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yi = — 1. Thus, we set y\ = — 1, y 2 = —2, and y% = —3. Satisfactory fits can be obtained 
after a cutoff for small system sizes L < 18, which yield fi c = 1.75682(2), y t = 1.51(1) ~ 
3/2 and Q = 0.823(2). It is interesting to observe that, without logarithmic corrections, 
satisfactory fits can include data for rather small sizes and the exponent y t = 1.51(1) agrees 
well with the exact value 3/2. This suggests that logarithmic corrections are small in the Q 
data, and thus the fitting results for Q are more or less reliable. 
We then fitted the x data by 

X = Xo + L 2y ^ 2 [a + ai (fi- n c )V* + a 2 (fi - p c )L 2yt 
+ b t L yi + b 2 L y2 + di/i - p c )L yt+y '} , (11) 

where Xo stems from the regular part of the free energy, which acts in Eq. (10) as a correction 
term with exponent y r = 2 — 2yh- The correction exponents were also set as y\ = — 1 
and y 2 = —2. After a cutoff for small systems L < 20, we obtain /i c = 1.75683(1) , 
y t = 1.489(9) w 3/2 and y h = 1.8748(8) ps y. The estimate of fx c is consistent with that 
from Q, and those for y t and yh agree with the exact values. If the exponent y t is fixed at 
3/2, one has /i c = 1.75683(1) and y h = 1.8743(7) after discarding the data for L < 18. 
The data for the particle density p were fitted by 

P = Po + Pi (A* - Pc) + L yt - d [a + ai (n - p c )L yt +a 2 (/i- p c )L 2yt 

+b 1 L' 1 + b 2 L- 2 + b 3 L' 3 } . (12) 

Satisfactory fits are obtained after a cutoff for small systems L < 16, and we have \i c = 
1.75682(2) , y t = 1.440(5) and p c = 0.180(4). The result /i c = 1.75680(3) agrees well 
with those obtained from magnetic quantities Q and x- However, the value y t = 1.440(5) 
significantly differs from 3/2. This might imply that, while additive logarithmic corrections 
are still small in energy-like quantities, multiplicative logarithmic corrections cannot be 
neglected. 

The finite-size scaling formula of the specific-heat-like quantities C e and C v reads 

C = r + n(i2- fi c ) + L 2yt - d [a + ai(/i - p c )L Vt + a 2 (p - p c )L 2yt + a 3 (p - Pc)L 3yt 

+ b 2 L~ 2 + b 3 L~ 3 + ciL 2 "- 1 (/1 - /i c )] . (13) 

In the actual fits, the terms r , arising from the regular part of free energy, cannot be 
distinguished from the correction term b\L 2yt ~ d ~ x = bi, and so is for r\ and c\. Thus, we 
simply set tq and r\ to be zero. The data for L > 18 are well described by Eq. (13). The 
fits for C v yield (i c = 1.75682(3) , y t = 1.468(7) and those for C e yield fi c = 1.75680(2) , 
y t = 1.470(5). Again, the estimates of /i c agree well with those from other quantities, while 
the values of y t differ significantly from 3/2. 

We have simulations at /j, — 1.756818, at criticality within the estimated error bars. Thus, 
we could analyze various quantities right at the critical point. The finite-size scaling behavior 
of x, P, and C e and C v at criticality is given by Eqs. (11), (12), and (13), respectively, by 
throwing out those /i-dependent terms. The estimates of the associated critical exponents 
from these simplified analyses are consistent with those from the aforementioned fits. Figures 
3 and 4 show the critical x an d C e data, respectively. 
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IV. DISCUSSION 



In the language of the lattice gas systems, we formulate a Swendsen- Wang-like version 
of the geometric cluster algorithm that has already found many applications [18,22]. Since 
our formulation is in line with the Swendsen- Wang algorithm for the ferromagnetic Potts 
model, we expect that it will help the reader to understand the geometric cluster method. 

We then study the triangular lattice gases with the first- and the second-neighbor ex- 
clusion, using a combination of the Metropolis and the geometric cluster algorithm. The 
estimated critical point ji c = 1.75682(2) significantly improves over the existing result 
1.7599; to our knowledge, no report has been published yet for the critical particle den- 
sity p c = 0.180(4). The excellent agreement between the exact values and the numerical 
estimates yt = 1.51(1) and y^ = 1.8743(7) give strong support for the expectation that the 
model is in the 4-state Potts universality class. On the other hand, the fitting results from 
energy-like quantities imply that, although additive logarithmic corrections might be small, 
multiplicative logarithmic corrections cannot be neglected at least in energy-like quantities. 

The fitting results are summarized in Table II. 
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FIGURES 




FIG. 1. The triangular lattice and its four sublattices. For chemical potential /j, — ► oo, one of 
the sublattices is fully occupied. 
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FIG. 2. Ratio Q versus (i for various system sizes. The straight line segments, simply connecting 
the data points, are for illustration purpose. 
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FIG. 3. Quantity x/L 2 at fi = 1.756818 versus L 2yh 4 , with y^ = 15/8. The statistical error 
bars are smaller than the size of the data points. The dashed line is just for illustration purpose. 




FIG. 4. Specific-heat-like quantity C e at \i = 1.756818 versus L 2yt ~ 2 , where the value of 
yt = 1.470 is taken form the fit. The statistical error bars of the data points are in the same 
order of their size. The dashed line is for illustration purpose. 
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TABLES 



TABLE I. Energies £\ and £2 of a building block in Eq. (2). The associated bond weight for 
K < in the geometric cluster algorithm is also given. In the "Example" , the upper two sites are 
i and i + 1, and the lower are — i and — i — 1. The filled (empty) circle represents the presence 
(absence) of a particle. 



Case 


1 


2 


3 


4 


5 


6 


7 




particle 


1 particle 




2 particles 


3 particles 


4 particles 


Example 


0-0 


•-0 


•—0 


•-• 


•—0 


•-• 


•-• 




0-0 


0-0 


•—0 


0-0 


0— • 


•—0 


•-• 


Si 











-K 





—K 


-2K 


£2 














-K 


-K 


-2K 


V 














e~ K -l 









TABLE II. Fitting results for various quantities. Symbol L m i n is the minimum system size for 
which the Monte Carlo data are included in the fit. 



Quantity L min fi c p c y t y h 



Q 


20 


1.75682(2) 




1.51(1) 




X 


20 


1.75683(1) 




1.489(9) 


1.8748(8) 


p 


18 


1.75680(3) 


0.180(4) 


1.440(5) 




£ 


12 


1.75680(2) 




1.48(2) 






18 


1.75680(2) 




1.470(5) 




C v 


18 


1.75682(3) 




1.468(7) 




Previous 




1.7599 




1.400 


1.885 
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